
if (free_childcare)
    P_c_min_Y_c=P_c.*min_Y_c;
    
    if ~no_response
        W=W0+P_c_min_Y_c;
    end
end

% unit cost of investment
if (Cobb_Douglas)
        
    % the share parameters already sum to one
    share_g=a_g;
    share_m=a_m;
    share_f=a_f;
    share_y=a_y;
    
    pbar=(p./share_g).^share_g.*(P_c./share_y).^share_y.*(W_m./share_m).^share_m.* ...
         (W_f./share_f).^share_f;
     
    pbar_single=(p./share_g).^share_g.*(P_c./share_y).^share_y.*(W_m./share_m).^share_m;
    pbar(married==0)=pbar_single(married==0);

    if (free_childcare)
        % expenditure shares out of home input expenditure
        share_g_h=a_g./(a_g+a_m+a_f);
        share_m_h=a_m./(a_g+a_m+a_f);
        share_f_h=a_f./(a_g+a_m+a_f);        
        
        
        pbar_h=(p./share_g_h).^share_g_h.*(W_m./share_m_h).^share_m_h.* ...
             (W_f./share_f_h).^share_f_h;
        
        pbar_h_single=(p./share_g_h).^share_g_h.*(W_m./share_m_h).^share_m_h;
        pbar_h(married==0)=pbar_h_single(married==0);
    
    end
    
else

    % tau_m/g
    Phi_m=(a_g./a_m.*W_m./p).^(1/(rho-1));
    
    % tau_f/g
    Phi_f=(a_g./a_f.*W_f./p).^(1/(rho-1));
    
    % define something for single and married
    aux=a_m.*Phi_m.^rho+a_f.*Phi_f.^rho+a_g;
    aux_single=a_m.*Phi_m.^rho+a_g;
    aux(married==0)=aux_single(married==0);
    
    % Y_c/g
    Phi_c=(a_g.*a_h./a_y.*P_c./p.*aux.^(gma/rho-1)).^(1/(gma-1));
    
    % g/X
    Theta_g=(a_h.*aux.^(gma/rho)+a_y.*Phi_c.^gma).^(-1/gma);
    
    pbar=Theta_g.*(p+P_c.*Phi_c+W_m.*Phi_m+W_f.*Phi_f);
    pbar_single=Theta_g.*(p+P_c.*Phi_c+W_m.*Phi_m);
    pbar(married==0)=pbar_single(married==0);
    
    if (free_childcare)
        % g/X_h
        Theta_g_h=aux.^(-1/rho);
        
        pbar_h=Theta_g_h.*(p+W_m.*Phi_m+W_f.*Phi_f);
        pbar_h_single=Theta_g_h.*(p+W_m.*Phi_m);
        pbar_h(married==0)=pbar_h_single(married==0);
        
        % expenditure shares out of home input expenditure
        share_m_h=Phi_m.*Theta_g_h.*W_m./pbar_h;
        share_f_h=Phi_f.*Theta_g_h.*W_f./pbar_h;
        share_g_h=Theta_g_h.*p./pbar_h;

        share_f_h(married==0)=0;
        
    end
    
    % expenditure shares
    share_m=Phi_m.*Theta_g.*W_m./pbar;
    share_f=Phi_f.*Theta_g.*W_f./pbar;
    share_g=Theta_g.*p./pbar;

    share_f(married==0)=0;
    share_y=1-share_m-share_f-share_g;

end

K=alpha.*beta.^(T-age+1).*delta2.^(T-age).*delta1;

% annual consumption, investment and leisure
sum_cshare=1+psi_m+psi_f+K; % sum of consumption share parameters

c=1./sum_cshare.*W./p;
X=K./sum_cshare.*W./pbar;
l_m=psi_m./sum_cshare.*W./W_m;
l_f=psi_f./sum_cshare.*W./W_f;

% per child weekly investment
X=X/52;

% investment expenditure (weekly per child)
pbar_X=pbar.*X;
p_g=share_g.*pbar_X;
W_m_tau_m=share_m.*pbar_X;
W_f_tau_f=share_f.*pbar_X;
P_c_Y_c=share_y.*pbar_X;

% investmenet quantity
g=p_g./p;
tau_m=W_m_tau_m./W_m;
tau_f=W_f_tau_f./W_f;
Y_c=P_c_Y_c./P_c;

% weekly hours worked
hrs_m=100-tau_m-l_m/52;
hrs_f=100-tau_f-l_f/52;

if free_childcare

    if no_response
        
        Y_c=Y_c+min_Y_c/52;
        
        if Cobb_Douglas
            X=tau_m.^a_m.*tau_f.^a_f.*g.^a_g.*Y_c.^a_y;
            X_single=tau_m.^a_m.*g.^a_g.*Y_c.^a_y;
        else
            X=((a_h.*(a_m.*tau_m.^rho+a_f.*tau_f.^rho+a_g.*g.^rho).^(gma/rho)+a_y.*Y_c.^gma).^(1/gma));
            X_single=((a_h.*(a_m.*tau_m.^rho+a_g.*g.^rho).^(gma/rho)+a_y.*Y_c.^gma).^(1/gma));
        end
        
        X(married==0)=X_single(married==0);
            
    else
        
        % check if Y_c>=min_Y_c
        for i=1:N

            min_x=1.d-3;
            if Cobb_Douglas
                max_x=min_Y_c(i)*P_c(i)/(pbar(i)*a_y(i));
            else
                max_x=min_Y_c(i)/(Phi_c(i)*Theta_g(i));
            end
            
            if (insample_baseline(i)==1 & edu_m(i)==1 & X(i)*52<max_x)
                
                x2=fzero(@foc_free_childcare,[min_x max_x]);
                
                % these are annual amounts
                X(i)=x2;

                if (Cobb_Douglas)
                    E=pbar_h(i)*(X(i)/(min_Y_c(i)^(a_y(i))))^(1/(1-a_y(i)));        
                else
                    E=pbar_h(i)*((X(i)^gma-a_y(i)*min_Y_c(i)^gma)/a_h(i))^(1/gma);
                end
                
                c(i)=(W0(i)-E)/(1+psi_m(i)+psi_f(i));
                l_m(i)=psi_m(i)*c(i)/W_m(i);
                l_f(i)=psi_f(i)*c(i)/W_f(i);
                
                % weekly amounts from now on
                X(i)=X(i)/52;
                E=E/52;
                
                % investment expenditure
                p_g(i)=share_g_h(i)*E;
                W_m_tau_m(i)=share_m_h(i)*E;
                W_f_tau_f(i)=share_f_h(i)*E;
                P_c_Y_c(i)=0;
                
                % investment quantity
                g(i)=p_g(i)/p(i);
                tau_m(i)=W_m_tau_m(i)/W_m(i);
                tau_f(i)=W_f_tau_f(i)/W_f(i);
                Y_c(i)=min_Y_c(i)/52;

                % hours worked
                hrs_m(i)=100-tau_m(i)-l_m(i)/52;
                hrs_f(i)=100-tau_f(i)-l_f(i)/52;
            end
        end
    end
end

% annual quasi difference in log achievement
qd_logPsi=delta1.*log(X.*52)+X_theta*phi_theta;

%% sample selection

if ((decomposition & ~baseline) | (price_reduction & ~quo) | free_childcare)
    insample=insample_baseline;
else
    insample=(insample_data==1 & hrs_m>0 & (married==0 | (married==1 & hrs_f>0)));     
end

qd_logPsi(insample==0)=nan;

p_g(insample==0)=nan;
W_m_tau_m(insample==0)=nan;
W_f_tau_f(insample==0)=nan;
P_c_Y_c(insample==0)=nan;

tau_m(insample==0)=nan;
tau_f(insample==0)=nan;
g(insample==0)=nan;
Y_c(insample==0)=nan;
X(insample==0)=nan;
pbar(insample==0)=nan;

%% model shares
share_m(insample==0)=nan;
share_f(insample==0)=nan;
share_g(insample==0)=nan;

for m=1:2
    for e=1:2
        
        eval(['share_model(m,e,1)=mean(share_m(ind' int2str(m) int2str(e) '),''omitnan'');']);
        if (m==2)
            eval(['share_model(m,e,2)=mean(share_f(ind' int2str(m) int2str(e) '),''omitnan'');']);
        end
        eval(['share_model(m,e,3)=mean(share_g(ind' int2str(m) int2str(e) '),''omitnan'');']);
    end
end

%% model levels
pbar_X(insample==0)=nan;
hrs_m(insample==0)=nan;
hrs_f(insample==0)=nan;

working_m(insample==0)=nan;
working_m(insample==1)=hrs_m(insample==1)>0;

for m=1:2
    for e=1:2

        eval(['level_model(m,e,1)=mean(tau_m(ind' int2str(m) int2str(e) '),''omitnan'');']);        
        eval(['level_model(m,e,2)=mean(hrs_m(ind' int2str(m) int2str(e) '),''omitnan'');']);
        if (m==2)
            eval(['level_model(m,e,3)=mean(hrs_f(ind' int2str(m) int2str(e) '),''omitnan'');']);
        end
        
    end
end

%% detailed outcomes
if (decomposition)
    
    outcome=zeros(14,2);
    doutcome=zeros(14,2);
    
    for e=1:2
        j=0;
        
        for m=1:2
            
            % expenditure
            j=j+1;
            eval(['outcome(j,e)=mean(pbar_X(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);

            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % price
            j=j+1;
            eval(['outcome(j,e)=mean(pbar(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);

            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end
            
            % total quantity            
            j=j+1;
            eval(['outcome(j,e)=mean(X(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % mother's time
            j=j+1;
            eval(['outcome(j,e)=mean(tau_m(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % mother's time expenditure
            j=j+1;
            eval(['outcome(j,e)=mean(W_m_tau_m(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % mother's time expenditure share
            j=j+1;
            eval(['outcome(j,e)=mean(share_m(ind' int2str(m) int2str(e) ...
                  '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)-outcome(j,1))*100;
            end
            
            % quasi difference in log achievement
            j=j+1;
            eval(['outcome(j,e)=mean(qd_logPsi(ind' int2str(m) int2str(e) ...
                  '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)-outcome(j,1))*100;
            end
        end
    end
    
    % variance decomposition
    vars=zeros(3,4);
    
    % put back data for expenditure
    pbar_X=pbar_X_data;
    pbar_X(insample==0)=nan;
    X=pbar_X./pbar;
    
    vars(1,1:3)=[var(log(pbar_X),'omitnan') var(log(pbar),'omitnan') var(log(X),'omitnan')];
    vars(2,1:3)=[var(log(pbar_X(married==0)),'omitnan') var(log(pbar(married==0)),'omitnan') var(log(X(married==0)),'omitnan')];
    vars(3,1:3)=[var(log(pbar_X(married==1)),'omitnan') var(log(pbar(married==1)),'omitnan') var(log(X(married==1)),'omitnan')];

    for m=1:3
        vars(m,4)=(vars(m,1)-vars(m,2)-vars(m,3))/(2*sqrt(vars(m,2)*vars(m,3)));
    end
elseif (free_childcare)

    outcome=zeros(20,2);
    doutcome=zeros(20,2);

    P_c_min_Y_c(insample==0)=nan;
    min_Y_c(insample==0)=nan;
    
    % calculate achievement at age 13
    logPsi=zeros(N,1);
    for t=6:T+1
        logPsi=delta2*logPsi+qd_logPsi;
    end
    
    for e=1:2
        j=0;
        
        for m=1:2
            
            % public childcare expenditure
            j=j+1;
            eval(['outcome(j,e)=mean(P_c_min_Y_c(ind' int2str(m) int2str(e) ...
                ')/52,''omitnan'');']);

            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % public childcare quantity
            j=j+1;
            eval(['outcome(j,e)=mean(min_Y_c(ind' int2str(m) int2str(e) ...
                ')/52,''omitnan'');']);

            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end
            
            % total quantity            
            j=j+1;
            eval(['outcome(j,e)=mean(X(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % mother's time
            j=j+1;
            eval(['outcome(j,e)=mean(tau_m(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % father's time
            j=j+1;
            eval(['outcome(j,e)=mean(tau_f(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % goods
            j=j+1;
            eval(['outcome(j,e)=mean(g(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % childcare
            j=j+1;
            eval(['outcome(j,e)=mean(Y_c(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % mother working
            j=j+1;
            eval(['outcome(j,e)=mean(working_m(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end
            
            % qd_logPsi
            j=j+1;
            eval(['outcome(j,e)=mean(qd_logPsi(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);

            % logPsi
            j=j+1;
            eval(['outcome(j,e)=mean(logPsi(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
        end
    end
end

